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Abstract 

With a view to a better understanding of the influence of atomic quantum delocalisation effects 
on the phase behaviour of water, path integral simulations have been undertaken for almost all of 
the known ice phases using the TIP4P/2005 model, in conjunction with the rigid rotor propagator 
proposed by Miiser and Berne [Phys. Rev. Lett. 77, 2638 (1996)]. The quantum contributions 
then being known, a new empirical model of water is developed (TIP4PQ/2005) which reproduces, 
to a good degree, a number of the physical properties of the ice phases, for example densities, 
structure and relative stabilities. 
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I. INTRODUCTION 



"Water, water, every where..." goes the poet Samuel Taylor Coleridge's The Rime of 
the Ancient Mariner, which provides a magnificent resume of our reason for studying this 
ubiquitous material. Many volumes have been written about water and ice (to cite just a 
few^"^^^), and a good deal more await writing, before we fully understand this enigmatic 
molecule. 

Currently the point has been reached where many properties, including the global phase 
diagram of water and the ice phases, can be reproduced qualitatively (and in some cases, 
quantitatively) using little more than a simple empirical model^. However, there are several 
aspects of water where our knowledge, and thus our understanding, is far from complete. 
One such aspect is the high pressure/temperature region of the phase diagram, where the 
precise location of the melting curves is still yet to be agreed upon due to the difficult 
nature of the experiments. For example, it is an open question as to whether water becomes 
super-ionic in this regiort^*^. In one of the ice polymorphs, ice X, the notion of a water 
molecule even becomes lost, the protons being shared equally between oxygen atoms^i^^. 
The low temperature region of the phase diagram is also extremely interesting, where a host 
of 'anomalous' or atypical trends are also present. Examples are the well known maximum 
in density at 3.984 Celsius, a minimum in the isothermal compressibility at 46.5 Celsius, and 
an unusual variation of the diffusion coefficient with pressure. These trends are especially 
apparent in super-cooled water where one can also find a minimum in both the density^ 
and a dynamic transition to Arrhenius behaviour for the diffusion coefficientp^ii^. It has 
been suggested that many of the anomalous properties of water at low temperatures could 
be understood by an hypothesised second critical poinir^^^ buried deep within "no-mans 
land"~, a region of the phase diagram inaccessible to experiment. If this is so, it would go 
a long way to explaining another feature of water; its capacity to form several amorphous 
phases (glasses) at low temperatures. 

In elucidating the origin of these anomalies, computer simulations have played a promi- 
nent role, for example their part in the proposal of a second critical point in water— 
using a simple empirical model. Classical computer simulations do, however, have their 
limitations. There are certain systems, water being one of them, where quantum effects are 
significantp^i22. As an example, let us examine the difference in temperature between the 
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melting point and the temperature of maximum density. For H2O this amounts to 3.984K, 
whereas for D2O it is 7.365K. From the point of view of the Born-Oppenheimer approxima- 
tion the potential energy surface (PES) is independent of the isotope considered. Thus the 
different behaviour of these isotopes is due to how the molecules react to this PES. This is 
known as an atomic quantum delocalisation effect. In this particular case the origin of the 
differences, both structural and dynamical, is in good part due to the quantum nature of 
the hydrogen protons and the strength of the hydrogen bond. Another example is the self- 
diffusion coefficient, which increases by more than 50% in a quantum system with respect 
to classical molecular dynamics simulationa^iiSI. 

The overall structure of water is that of an asymmetric top, which is to say that all 
three principal moments of inertia are distinct. What is particularly interesting is that 
since hydrogen is the lightest atom, the rotational moments of inertia are small enough 
to show marked quantum behaviour. Thus water has significant quantum effects even at 
room temperature. The importance of these quantum effects increases as the temperature 
is lowered. For the ice phases these effects are expected to be significant, especially at 77K 
where many experiments on ice are frequently performed using liquid nitrogen. Thus far 
there has been relatively little work on these effects for ice, and almost all of the work that 
has been published has focused on ice I;^— The objective of this publication is to 
quantify the size of these effects in all of the ice phases, apart from that of ice X, which 
cannot be described by the rigid models used in this work. 

These atomic quantum delocalisation effects will be studied using the empirical 
TIP4P/2005 model^. Over the last few years a number of the present authors have un- 
dertaken extensive simulation studies examining the performance of a number of commonly 
used models for water, in particular the TIP3P, T1P4P, T1P5P and SPC/E models^. The 
principal findings have been that the TIP3P— , TIP5P— and SPC/E^° models experience 
difficulties when it comes to describing the global phase diagram of water and the ice phases. 
However, the TIP4P model does indeed provide a qualitatively correct phase diagram. Based 
on this finding, the TIP4P model was re-parameterised in order to improve the quantitative 
representation, leading to the T1P4P/2005 model^. It has since been found that this model 
also provides a good description of the maximum in density of liquid water and its variation 
with pressure^, of the compressibility minima^^, the surface tensionS^ the vapour liquid 
equilibria^, the critical properties^, the equation of state at high pressures^^, the diffusion 
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coefficienlj^, and the viscosity^. 

That said, the model was parameterised for classical simulations, so the introduction 
of atomic quantum delocalisation effects, although improving the qualitative description, 
will cause a deterioration in the quantitative description. In the first stage of this research 
we shall analyse the impact of atomic quantum delocalisation effects on the properties of 
the ice phases using this potential. That will elucidate where, and how, atomic quantum 
delocalisation effects modify the properties of water with respect to the classical limit. These 
differences then known, we provide a re-parameterised version of the TIP4P/2005 model 
which we shall call the TIP4PQ/2005 model, the Q indicating that this model is suitable 
for quantum simulations. As was pointed out by Morse and Rice^^ as well as by Whalley 
"...effective potentials that are used to simulate water ought to be tested on the many phases 
of ice before being treated as serious representations of liquid water"—. 



II. METHODOLOGY 



Simulations were performed using the path integral formulation, which permits us to 
study the quantum effects related to the finite mass of the atoms (in many quantum chem- 
istry calculations, the electrons are treated as being quantum, however the nuclei are treated 
as classical point masses). A particularly elegant technique for studying quantum effects in 
many body systems is that of path- integral Monte Carlo (PIMC). There are many good in- 
troductions concerning PIMC in the literature^i^»^'^i^, here we shall focus on the aspects 
most pertinent to the simulations we have performed. 

Water is, of course, a flexible molecule. For path integral simulations one generally 
requires the number of Trotter slices, P, to be^ 



where cUmax is the 'fastest' frequency present in the system in question. In water the in- 
tramolecular vibrations are of the order of co'max/27rc ~ 4000 cm~^ which leads to P > 20. 
Using the rigid body approximation for water the fastest motion now becomes the libration, 
with a frequency of < 900 cm^^, thus reducing P to around 5-6. This represents a substan- 
tial reduction in the computational overhead associated with traditional PIMC calculations 
(although new techniques have recently been developed by Manolopoulos et al. to increase 
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the efficiency of flexible molecule PIMG^). It must be said that by choosing to use a rigid 
model, one precludes the ability to study some aspects of water, such as the high frequency 
region of the infra-red adsorption spectrum^i^. The infra-red spectrum of water and ice 
can be divided up into two distinct regions. Above ~ 900 cm^^ one has the contribution 
associated with the intramolecular degrees of freedom of bending and stretching. Below 
^ 900 cm~^, as previously mentioned, one has the section that corresponds to translational 
and librational movements, and are mostly due to inter-molecular forces. Quantum contri- 
butions to the Helmholtz energy (A) within a perturbative treatment for a rigid asymmetric 
top are given by^: 

cyclic ^ ' 

A good proportion of the quantum effects in water are due to the strength of the hydrogen 
bond, along with a particularly small inertia tensor. It is this that lends importance to 
the torque (F) terms found in the above equation, which results in the appearance of the 
librational band. In contrast, this region for a molecule such as SO2, where no such hydrogen 
bonding is present, is far less important. By using the path integral formulation for a rigid 
model we shall be studying atomic quantum delocalisation effects in the influential region 
encountered below ^ 900 cm~^. In a study of the phonon density of states for ice Dong 
and Li^^ showed that the rigid TIP4P model does a reasonable job of reproducing this 
low frequency section of the spectrum. Even given the fact that intramolecular effects are 
important, it is surely the case that a rigid body path- integral study is more physically 
realistic than a purely classical study, which neglects all atomic quantum delocalisation 
effects. Such an approach has been adopted in a number of studies, using for example 
the SPC/E model2^. In view of this, and given the success that the TIP4P/2005 model 
has had in describing the ice phases classically, the rigid TIP4P/2005 model is the natural 
candidate for a preliminary study of atomic quantum delocalisation effects in ices. Given 
that the TIP4P/2005 model is a rigid asymmetric top, we shall first present the path integral 
description of a rigid rotor. 
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A. Path integrals for a rigid molecule 



The coordinates used to describe a rigid molecule are rifii, where ri represents the centre 
of mass and fli = {4>i,0i,Xi) represents the Euler angles that fix the molecule orientation. 
The Hamiltonian of a rigid asymmetric rotor can be written in the forn>^: 



= T"^ + T 



-irot 



u. 



(3) 



where T^^'^ represents the kinetic energy operator associated to the centre of mass translation, 
U appears as a potential energy operator that is a function of the coordinates rifii, and the 
rotational kinetic energy operator is given by^: 

3 f2 



T 



rot 



where Li are the components of the angular momentum operator and Jj are the moments 
of inertia of the molecule referred to its fixed body frame. We assume, without loss of 
generality, that the moment of inertia tensor is diagonal in the chosen fixed body frame. 

In the path integral formulation, the partition function, Qi, of a rigid molecule may be 
expressed by a factorisation of the density matrix into P factors, so that each quantum 
particle is described by a ring of P replicas or 'beads', 

p p 



QM 



lim 

P-*oo 



t.t+1 

Pi 



(5) 



i=l 



t=l 



where j3 
by^: 



l/ksT is the inverse temperature, and the propagator pY^^ is approximated 



pY^\P/p) 



r\n\\exp -PU/2P 



exp 



-(3{T''^ + T'°')/P 



exp 



-PU/2P 



(6) 



The propagator satisfies the cyclic condition that bead P + 1 corresponds to bead 1. This 
rigid molecule propagator is built up of three factors, a potential energy component, a 
translational component, and a rotational component: 



t,t+i, 
Pi ' 



L/o/p\^ t,t+l t,t+l t,t+l 
\P I r) ~ Ppot,lPtra,l Prot,l • 

The potential energy component is given by^ 

exp 



t,t+i 

Ppot,l 



2P ^ 



(7) 



(8) 
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(9) 



where = U{r\n\) is the potential energy of the rephca t of the molecule. The translational 
component is given by^ 

3/2 

where M is the total mass of the rigid molecule. The two previous equations are well known 
and are commonly used as the so-called primitive approximation in path integral studies of 
simple fluids. The rotational propagator between t and t + 1 is given by^: 



,-N(r;|exp(-,f-/P)|rS«)^(i|3)'^xp 



pSi = {^il'^^p [-i3T'°'/p\ |n;«^ . (10) 

In an important piece of work Miiser and Berne^i^ have shown that the rotational contri- 
bution to the propagator between the replicas t and t + 1 of a rigid molecule i is exactly 
given by 

oo J J / fDpJ]\l\ 

J=OM=-Jk=-J V / 

where 

/• rii^ = ^^mm(^T"') cos[M(0r^ + Xr^)]l<fl' (12) 
Here dj^M^i'*^^) Wigner functions and l/l^^-*! are the coefficients of the expansion of 
the eigenstates of the asymmetric top in a basis formed by the eigenstates of the symmetric 
top. E^^^^^ are the eigenvalues of the energy of the asymmetric top. The quantum numbers 
J and M provide the values of the total angular momenta of the asymmetric top and the 
value of its z component in the laboratory frame. The number K is not a true quantum 
number, in the sense that it does not provide the value of any physical observable, but rather 
is an index used to label the (2 J + 1) energy levels that are obtained for each value of J. 
The angles 6'*'*'*'^, 0*'*"^^ and Xi'*^^ are the Euler angles of the replica t + 1 of molecule i 
expressed in the body frame fixed in the replica t of the same molecule i. Note that the 
rotational propagator depends solely on two variables, and 0*'*^^ + Xi'*^^- Obviously 

to determine the value of the rotational propagator one must first determine the (2J + 1) 
energy levels of the asymmetric top for each value of J. This can be obtained from the 

43 . The coefficients lAl-^f^"*! are 

' I KM I 



(2 J + 1) eigen- values, E^^^^\ of the matrix given in Ref. 
the eigen-vectors associated with these eigen-values. It is computationally convenient to 
calculate the rotational propagator plot,i {Of^~^^ , + Xi'*"'^^) for a gi'id of values of the 



angles 6'*'*'*'^ and <^*'*^^ + xV~^^ for sach value of /3/P to be used, and save this data prior to 
the simulations. The value of the rotational propagator for any given i?*'*"*"^ and 0*'*^^ + 
can then be estimated using a linear interpolation algorithm from this tabulated data. 

B. Path integrals for an ensemble of rigid molecules 

Once the translational and rotational propagators are known for a rigid molecule one can 
calculate the partition function for a set of interacting molecules. Let us assume that we 
shall be using a pair-wise potential Uij such that the potential energy of the replica t of the 
system is 

u' = j2J2''^Myp^i^',)- (13) 

i j>i 

Now the canonical partition function, Qn, of an ensemble of molecules described with P 
beads is given by: 

1 / MP \3JVP/2 . ^ N P 

-p (-H E E - ^rr - ^ E n n ■ 

V ' i=l t=l t=l / i=l t=l 

As can be seen in Eqs. f[T^ and ([HD, each replica t of molecule i interacts: (a) with 
the molecules that have the same index t via the intermolecular potential Uij; (h) with 
replicas t — 1 and t + 1 of the same molecule i via a harmonic potential whose coupling 
parameter depends on the mass of the molecules, M, and on the inverse temperature /?; and 
(c) with replicas t — 1 and t + 1 of the same molecule through the terms p^^l^l and p*ot^/ 
which incorporate the quantisation of the rotation, which in turn depends on the relative 
orientation of replica t with respect to t — 1, and t + 1 with respect to t. 
Let us define an energy U' as: 

MP ^ ^ 1 ^ 

and the total orientational propagator Pj-ot as: 

N P 

Pro. = \{\{pYot} (16) 

i=l t=l 
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Within a Monte Carlo simulation one generates a new configuration starting from a previous 
configuration. The probability of accepting this new configuration, Paccept, is given by 



Paccept = mm 



pnew 

l,exp(-/3(f/L^-f/:j)-g;^ 

-'rot J 



(17) 



It is worthwhile making two observations about the orientational propagator between a pair 
of contiguous beads pH^l ■ Firstly, it must be positive in order to be used in the Metropolis 
acceptance criteria, which is indeed the case. Secondly, the maximum in the orientational 
propagator is achieved when 9 = and 4>+x = 0. It is found that at high enough temperature 
the propagator decays to zero relatively quickly as the values of 6 and (p + x increase. The 
orientational propagator can also be expressed as an auxiliary energy by defining Ui^aux such 
that 



u. 



t,t+i _ 1 , t,t+i 



InpJ:,?/ (18) 



t^aux p 

Ui^aux has a minimum at = and + x = and increases quickly as a function of the 
variables 6 and </> + X- Pmt can now be written as 

Prot = exp(-/5f/,,J = exp ( -/3 uX£\ . (19) 

\ i=i t=i / 

Using this auxiliary energy the Metropolis criteria can be now written as : 
Paccept = nim 

This expression helps us to clarify the role of the orientational propagator; it can be viewed 
as a potential that forces two contiguous beads, t and to adopt similar orientations (this 
corresponds to the minimum of the auxiliary potential) with an energetic penalty when they 
adopt different orientations. This is analogous to the role played by the harmonic springs 
connecting the centre of masses of the molecules in Eq. f|T5l) . 
The internal energy can now be calculated from: 

i? = - — (21) 

Qn 013 ^ > 

It can be shown that substituting the value of the canonical partition function in this ex- 
pression results in 

E = Ktra + i^rot + t/, (22) 
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where : 

N P 



i=l t=l 

1 sr-^sr-^ l^J=0 l^M=-J l^K=-J J J M k k ^ 



'p^k 



1=1 t=l Prot,i 

U=i^j:U^y (23) 

As with the rotational propagator, the numerator of Krot in Eq- l23] was calculated prior to 
the simulations for a grid of the variables 6 and (p + x ^"^^ subsequently saved in tabular 
form. 

When performing simulations of solids it is more convenient to perform the simulations 
in the NpT ensemble. The partition function for the NpT ensemble can be calculated using: 

Qnpt = aJ dV exp{-(3pV)QN (24) 

where A is a constant with units of inverse volume that makes Q dimensionless. Its value 
affects the Helmholtz energy function, but not the configurational properties. 



C. Simulation details 



In this work path integral Monte Carlo simulations are undertaken for the TIP4P/2005 
model for fourteen of the fifteen known ice phases. One of the most important variables when 
it comes to path integral simulations is the number of Trotter slices, or beads, (P) employed. 
If P = 1 then the simulation is classical. As P ^ oo then the quantum simulation becomes 
exact. Given the isomorphism between Trotter slices and the number of component 'beads' 
in a ring polymer—, one can easily see that the time required for a simulation scales with the 
number of Trotter slices used. For flexible models of water at 300K a typical number of slices 
is about P = 24^1^1^. However, if a rigid model is employed, the number of Trotter slices 
required can be reduced by about a factor of five^^*^. Previous studies for a rigid model 
of water at 300K found that a value of P = 5 provides good convergence^^i^. Thus in this 
work the number of Trotter slices times the temperature was maintained at PT ^ 1500. For 
the lowest considered temperature (77K) this corresponds to 20 beads. When computing 
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the asymmetric top eigen-energies and eigenvectors of water the OH distance and the H- 
0-H bond angle of the TIP4P/2005 model were used, which corresponds to the gas phase 
geometry of real water. The principal moments of inertia are computed using this geometry 
along with the masses of the hydrogen and oxygen atoms. Although the model has the 
negative charge on the site M, this site is massless and therefore it is only used to compute 
the potential energy of the system. 

In this work two models of water are studied, the TIP4P/2005 modeP^ and a re- 
parameterisation, which we shall call the TIP4PQ/2005 model, to 'compensate' for quantum 
effects. The parameters for both of these models are given in Table [11 The only difference be- 
tween these models is an increase in the charges on the hydrogen sites by 0.02e, along with a 
corresponding increase in the charge on the oxygen site. For both models the Lennard- Jones 
potential was truncated at 8.5A and long-range corrections were included. The TIP4P/2005 
model has been designed to be used with Ewald summations^i^ which is a well known 
technique to treat the long range electrostatic interactions. Ewald summation is more ap- 
propriate than the reaction field method when it comes to the simulation of solid phases. 
The real part of the Coulombic potential was truncated at 8.5A. 

The configurational space of the quantum system was sampled using a Monte Carlo code 
with four distinct types of trial moves: the displacement of a single bead of one molecule, 
rotation of a single bead of one molecule, translation of a whole ring, and rotation of all 
of the replicas of one molecule. A Monte Carlo cycle is defined as Monte Carlo moves, 
where the probability of attempting a translation or a rotation of a single bead is 30% each 
and the probability of attempting a translation of a whole ring or rotating all the replicas 
of a ring is 20% each. The maximum displacement or rotation in each type of movement 
was adjusted to obtain a 40% acceptance probability. When simulations were performed in 
the NpT ensemble, besides the A^ particle trial moves, one Monte Carlo cycle also includes 
an attempt to change the volume of the simulation box. The maximum volume change was 
adjusted so as to obtain a 30% acceptance probability. In general the simulations consisted 
of 30,000 Monte Carlo equilibration cycles, followed by a further 100,000 cycles for the 
accumulation of run averages. The number of molecules used in each of the phases are given 
in Table [TTl For the proton disordered ice phases the positions of the hydrogen atoms were 
generated in such a way as to produce a system that satisfies the so-called Bernal- Fowler ice 
rules^"^, and whose dipole moment as close as possible to zero. This was achieved using 
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the algorithm proposed by Buch et a/.— 

As mentioned, all simulations were performed in the isothermal-isobaric (NpT) ensemble. 
The implementation of the NpT ensemble in PIMC has already been discussed in previous 
works^i^. It is important to note that the Monte Carlo volume moves should be performed 
anisotropically, in order to allow the simulation box to 'relax' and obtain the true equilibrium 
unit cell of the model under consideration. In other words, the pressure on the simulation 
box should be hydrostatic; the pressure tensor is diagonal and each of the elements along 
the diagonal have the same value. If this is not the case the system will suffer stresses and 
the structure and thermodynamic properties will not reach their equilibrium values. This 
is achieved using the technique proposed by Parrinello and Rahman^i^!^ and extended to 
Monte Carlo by Yashonath and Rao^. Briefly, the shape of the simulation box is defined 
by a so-called iJ-matrix representing the Cartesian coordinates of the vectors defining the 
simulation box. Each of the individual components of the if-matrix are adjusted randomly, 
leading to changes in both the simulation box lengths and in the geometry. 

As a preliminary check that the Miiser and Berne propagator was implemented correctly 
the rotational energies were calculated for an isolated H2O molecule. In Fig. [T]the rotational 
energies computed from the exact expression of the quantum partition function of an asym- 
metric topS^ (with the appropriate rotational constants) are compared to those obtained 
from PIMC simulations. As can be seen the agreement is excellent. It should be noted that 
the present calculations do not include exchange effects. However, these are only expected 
to be relevant at temperatures below those that we have studied in this work. 

III. RESULTS 

A single state point has been simulated for each of the solid phases of water with the 
exception of ice X, which cannot be described by a rigid model^'^''. The results of these 
simulations are presented in Table [Tll By comparing the densities obtained from classi- 
cal TIP4P/2005 simulations to path integral simulations of the TIP4P/2005 model, which 
henceforth we shall denote as TIP4P/2005(pi), it is clear that the introduction of atomic 
quantum delocalisation effects reduces the density of the solid phase by about 0.02 g/cm^ 
for temperatures above 200K, and by ^ 0.03 — 0.04g/cm^ for temperatures in the range 
75-170K. Not surprisingly, quantum effects become increasingly evident as the temperature 
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is reduced. The various contributions to the total energy, E, are also tabulated. As far 
as the translational kinetic energy component, -ft'transiationai, is concerned one can observe 
an increase of about 10% for TIP4P/2005(pi) with respect to TIP4P/2005 (i.e. (3/2)i?T) 
at temperatures above 225K. As the temperature is lowered, this difference becomes 100%. 
This is approximately true for all of the ices. From these results one can conclude that the 
translational contribution to the heat capacity in quantum simulations is significantly less 
than the corresponding contribution in classical simulations. If one looks at the rotational 
kinetic energy contribution, -ft"rotationah the differences are exaggerated even further; ranging 
from about 100% for 'high' temperature ices, and increasing to 600% at low temperatures. 
From this it is clear that the quantum contributions are manifestly rotational in their nature, 
whilst translational effects are secondary in the solid phase. Within a perturbative treat- 
ment the quantum contribution to the Helmholtz energy function is related to the average 
of the forces divided by the masses for the translational contribution, and to the average of 
the torques divided by the principal moments of inertia for the orientational contribution^. 
The mass of water is almost the same as that of neon, however, the quantum effects are 
far more pronounced in water for the temperature range considered in this work—. The 
overwhelming reason for this difference is the strength and directionality of the hydrogen 
bond. This, as well as the fact that the moments of the inertia tensor are quite small due 
to hydrogen having a very low mass. The temperature dependence of the kinetic rotational 
energy is rather weak, so its contribution to the heat capacity is expected to be small. On 
the other hand the quantum contributions to the potential energy are of the order of 1 
kcal/mol at high temperatures, which increases to 1.5 kcal/mol at low temperatures. Thus 
there is a significant difference in E between the TIP4P/2005 and the TIP4P/2005(pi) re- 
sults, amounting to about 3 kcal/mol at low temperatures; half of which being due potential 
energy, and the other half kinetic. 

We shall now turn to the radial distribution functions. These histograms provide insights 
into the structure of a fluid on a molecular scale^"^. One of the first simulation studies of 
such functions for water using path integral simulations was undertaken by Kuharsky, Rossky 
and co-workers^i^i^i^ for the ST2 model. Given the low scattering factor of hydrogen, the 
oxygen-oxygen (goo) is the distribution function most accessible experimentally. Here we 
present the oxygen-oxygen radial distribution function for ices Ih, II and VI (Figures [2lll]) for 
classical TIP4P/2005 and TIP4P/2005(pi). For ice Ih the experimental radial distribution 
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function has also been plotted, using the data provided by Soper— at 220K. To the best of 
our knowledge as yet there are no experimental radial distribution functions available in the 
literature for ices II and VI. In Table IIIII details are given for specific points located along 
the oxygen-oxygen radial distribution function curves for ice Ih- On going from classical 
simulations to path integral simulations the location of the first two peaks shifts to slightly 
larger distances. Furthermore, there is a notable reduction in the height of these peaks when 
quantum contributions are incorporated. Similar findings have been published previously 
for water and for simulations of TIP4P(pi) of ice I^ by Hernandez de la Pena et al.— . This 
softening of the distribution functions goes hand-in-hand with the reduction in the density 
of the ices in the PIMC calculations. It is interesting to speculate whether the addition of 
the small (and somewhat unusual) first peak in the ice 1^ experimental data with the much 
larger second peak would place the simulation results in a more favourable light. 

A consequence of the third law of thermodynamics is that the coefficient of thermal 
expansion, a, tends to zero when the temperature goes to zero. Experimentally one finds 
that there is very little variation in the density of ice Ih in the temperature range 0-125K. 
Classical simulations are unable to capture this, as can be seen in the low temperature 



equations of state published in Ref 



74J . where the density of ice continues to increase as 



the temperature is lowered. Here we have performed simulations of TIP4P/2005(pi) for 
temperatures in the range 77-200K along the atmospheric pressure isobar for a number of 
ices. These results are presented in Table HVl In particular, the equation of state of ice I^ is 
plotted in Fig. [5] along with classical^ and experimental results^. One can see a dramatic 
reduction in the density between classical TIP4P/2005 and TIP4P/2005(pi) simulations. 
However, the most important difference is that the density is almost independent of the 
temperature below ^ 125K, in other words, a tends to zero. Given the fact that the 
TIP4P/2005 model was parameterised for classical simulations, it is no surprise that the 
TIP4P/2005(pi) results show a significant deviation from the experimental values. That 
said, the TIP4P/2005(pi) curve is, more or less, parallel to the experimental curve, strongly 
suggesting that a re-parameterisation of the TIP4P/2005 model could improve these results 
by shifting the TIP4P/2005(pi) curve to higher densities. It is worth mentioning that the 
lOOK state point for the TIP4P/2005(pi) model seems to be slightly more dense than the 77K 
state point. It has been suggested that there is a temperature of maximum density in the 
ice phase^"^, however, longer and more detailed simulations would have to be undertaken 
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to establish whether our results do indeed reflect this or not, given that this curvature could 
well be due to the statistical uncertainties in the simulation results. 

In 1984 Whalley estimated the thermodynamic properties of ices at OK. This estimate was 
made after analysing the experimental coexistence curves between ices at low temperatures^ 
and realising that at OK phase transitions occur with zero enthalpy change. By assuming 
that the volume and internal energy difference between ices is largely unaffected by pressure 
(a quite reasonable approximation) Whalley was able to estimate the energies and densities 
of ices at OK and zero pressure. Such a calculation is useful as it allows one to obtain an idea 
of the form of the phase diagram at low temperatures by examining the relative stability of 
the ice phases. Thus one can estimate the coexistence pressure between two ice phases at 
zero kelvin using the approximation 

(25) 

p=0 
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More recently a similar analysis was undertaken^ for a number of popular empirical models 
of water. For the SPC/E and TIP5P models ice II was found to be more stable than ice I/i, 
however, for TIP4P/2005 ice I/j, as is the experimental situation, was more stable than ice 
II. Here simulations were performed at 125, 100 and 77K for TIP4P/2005(pi) (for technical 
reasons PIMC simulations at OK are infeasible, given the number of beads required). As- 
suming that the heat capacity, Cp, follows the Debye law, i.e Cp oc T'^, then it follows that 
the enthalpy should scale as T'^. Note that the internal energy and enthalpies are almost 
indistinguishable at room pressure; the pV term is negligible compared to the internal en- 
ergy term. In Fig. [6] the internal energies from Table [IV] are plotted as a function of the 
temperature for TIP4P/2005(pi) and the estimated values at OK, obtained from a fit of the 
form E = a + bT^, are given in Table El The relative energies between ices obtained at 
OK from the extrapolation procedure described above are quite similar to those obtained 
from the simulations results at 77K. The inclusion of quantum effects consistently increases 
the energy at OK of the ice phases by ~ 3.5 kcal/mol. However, for ices II, HI, V and VI 
the relative energy remains largely unchanged; differing by only 0.1 kcal/mol from the 
classical values. The zero point energies of ices II, HI, V and VI are quite similar and are 
expected to have very little effect on the relative stability of the ice phases. This is not 
the case for ice I^, atomic quantum delocalisation effects destabilise ice Ih with respect to 
ice II, the difference now being ^ 0.26 kcal/mol. For example, for TIP4P/2005(pi) ice II 
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replaces Ih as the most stable ice phase at low temperatures. Given the fact that quan- 
tum effects stabilise ice II with respect to ice Ih implies that for the TIP3P, SPC/E, and 
TIP5P models the inclusion of atomic quantum delocalisation effects would further deterio- 
rate their phase diagrams; the ice Ih phase being stable only for large negative pressures and 
ice II dominating the low temperate atmospheric pressure isobar. An interesting question 
is precisely why ice Ih is more affected than the rest of the ices by these atomic quantum 
delocalisation effects. As discussed previously, within a perturbative treatment the effect of 
atomic quantum delocalisation effects can be expressed as the average of forces and torques 
on the molecules divided by their masses or principal moments of inertia. Since the mass 
and inertia tensors are the same, regardless of the ice phase considered, differences between 
ices must be due to differences in forces and torques between molecules. In all the ices each 
water molecule forms four hydrogen bonds with its nearest neighbours. For ice Ih, the four 
nearest neighbours form an almost perfect tetrahedron. However, for ices II, III, V and VI, 
the four nearest bonds form a distorted tetrahedron^, resulting in weaker hydrogen bonds 
(even though they are more dense than ice Ih)- It is the strength of the Ih hydrogen bonding 
that is showing up in the quantum contributions. 

The results presented thus far have elucidated the quantum contributions to the properties 
of the solid phases of water. The TIP4P/2005 model used in this study was originally 
parameterised to reproduce as faithfully as possible the experimental results for water using 
classical simulations. Thus in some implicit way, quantum contributions form part of the 
make-up of this model. It is no surprise that an explicit introduction of quantum effects will 
degrade the qualitative aspects of this model, which is exactly what we have seen in this work 
using TIP4P/2005(pi). We have witnessed that quantum effects decrease both the structure 
and the density of the ices as the temperature is lowered, and that they modify the relative 
stability of ices Ih and II. Originally the TIP4P/2005 model was created by examining the 
derivatives of the parameters of the model for a number of properties, and then, via a least 
squares fit, the optimum values for the parameters are obtained. These properties include 
the density and the coexistence lines obtained from values of the Helmholtz energy function. 
However, here we do not yet have access to the coexistence lines for the TIP4P/2005(pi) 
model so in developing the new TIP4PQ/2005 model a modest, and quite probably sub- 
optimal, change in the parameters was called for. 

There is a veritable plethora of classical empirical models for water in the literature. In 
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contrast, there is a paucity of quantum empirical models. It is worth making a mention of 
three of these quantum models; a re-parameterisation of a flexible version of the SPC/Fw 
model^, the second is a re-parameterisation of the rigid TIP5P model^, and the third is a 
series of flexible and polarisable potential models named TTM2-F— and TTM3-F— , obtained 
from fits to the potential energies of water clusters obtained from first principle calculations. 
For both the SPC and the TIP5P re-parameterisations the essential difference was that the 
dipole moment of the molecule was increased, whilst maintaining the remaining parameters 
of the potential constant. The basic idea is that since atomic quantum delocalisation effects 
reduce the density and internal energy of the system, increasing the charge is a simple way of 
're-compensating' for these changes, coaxing the model back to being its former self. It was 
with this in mind that the TIP4PQ/2005 model was created. The only difference between 
the TIP4P/2005 and the TIP4PQ/2005 models is in the dipole moment (see Table H]), which 
was increased from 2.305D to 2.38D. This was achieved by a 0.02e increase in the charge 
of the protons. Similar increases in the dipole moments of water (of about 0.08-O.lOD) 
were used in the aforementioned quantum versions of SPG^ and TIP5P models^. Such an 
increase in the charge may not be necessary in a flexible model where, as stated by Mahoney 
and Jorgensen, "...although quantum effects make the density behaviour of the rigid model 
worse, they improve the density behaviour of the flexible model."—. This interplay between 
an increase in the dipole moment and flexibility has also been commented upon by other 
authors^i^. Obviously this new model is only suitable for quantum simulations of water. 

In Table |Vl] the state points for the ice phases are recalculated using this new 
TIP4PQ/2005 model. When compared to the experimental ys,\^e^^^^^^^^^,93^4^ 
the results are really quite good over the whole range of temperatures and pressures. The 
average quadratic deviation between experimental and predicted densities (excluding ice 
VII) is 0.01 g/cm^ for the classical TIP4P/2005 model, which becomes 0.03 g/cm^ for the 
TIP4P/2005(pi) model. For the re-parameterised TIP4PQ/2005 model the quadratic de- 
viation is once again 0.01 g/cm^, recovering the situation for the classical model for the 
state points considered. In Table lylTl the unit cell parameters for the TIP4PQ/2005 model 
for a selection of ice phases have been provided and are also seen to be rather good when 
compared to the experimental values. 

In Fig [5] the equation of state for ice Ih is plotted. The TIP4PQ/2005 state points are 
equidistant from those of TIP4P/2005(pi), but they are now much closer to the experimental 
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values, with a deviation of around 0.005 g/cm^, which amounts to a difference of only 0.8% 
with respect to the experimental value. Given the curvature of the equation of state, in line 
with the third law of thermodynamics, and the small difference between the TIP4PQ/2005 
densities and the experimental results, leads us to believe that this is one of the best theoret- 
ical descriptions of ice Ih thus far seen in the literature. This is not to say that in the future 
this cannot be improved upon, for example via the inclusion of flexibility, polarisability etc. 
in the molecular model. In Fig. [7] the oxygen-oxygen radial distribution function of ice 1^ 
at 77K is compared to the experimental results of Narten^^, and the results are acceptable 
almost all the way up to 9A. The most notable difference can be seen in the height of the 
first peak; which drops from 9.37 for classical TIP4P/2005, down to 6.21 for TIP4PQ/2005, 
compared to 5.95 experimentally^^. 

In an analogous study to that for OK for TIP4P/2005(pi) the relative stability of ices I^,, 
II, III, V and VI at low temperatures has been tabulated in Tables IVl and \VTU\ and plotted 
in Fig. [HI As can be seen the relative energy between ice II and the remainder of the ices is 
similar to that of TIP4P/2005(pi). The most significant result is that for TIP4PQ/2005 ice 1^ 
regains its rightful place as the most stable ice phase. Experimentally the energy difference 
between 1^ and II is 0.014 kcal/mol, which for TIP4PQ/2005 becomes 0.04 kcal/mol. In 
Table [IXl results for the OK coexistence pressures, calculated using equation[25l are presented. 
It can be seen that both the energies (Table |V]) and the coexistence pressures (Table IIX|) for 
various transitions are substantially better than the values provided by classical simulations 
of the TIP4P/2005 model, in particular for the Ih-ll transition. This gives us confidence 
that the TIP4PQ/2005 could well produce a respectable global phase diagram in the future. 

IV. CONCLUSIONS 

This work addresses a series of physical properties of water that vary with the inclusion 
of atomic quantum delocalisation effects, which were introduced to the TIP4P/2005 model 
using path integral Monte Carlo simulations. Quantum simulations have been undertaken 
for all of the ice phases of water, with the exception of ice X, for the TIP4P/2005 model, 
and for the new TIP4PQ/2005 model. Using the Miiser and Berne propagator for rigid 
asymmetric tops, various properties of these ices have been examined. 

It has been found that the radial distribution functions become more 'washed-out' when 
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quantum effects are taken into account. In other words, the peaks become lower and wider 
and shift to shghtly larger distances. This goes hand-in-hand with a reduction in density for 
the quantum solid; by ~ 0.02 g/cm^ for temperatures above 150K, and ~ 0.04 g/cm^ below 
lOOK. 

If a classical empirical model is tailored to reproduce the experimental ice densities at a 
temperature close to the melting point, as the temperature is reduced the model will start 
to fail (such is the case, for example, of the TIP4P/2005 modet^). This is due to the fact 
that classical simulations are unable to satisfy one of the consequences of the third law of 
thermodynamics, namely that the coefficient of thermal expansion, a, tends to zero as the 
temperatures approaches zero Kelvin. It can be seen that the PIMC simulations now, to a 
good degree, correctly describe the low temperature physics of this model. 

The translational component of the kinetic energy bears a passing resemblance to the 
classical value of (3/2)i?T, whereas the rotational component is markedly larger. 

There is a particularly pronounced effect in the relative stabilities of ices Ih and II, where 
the stability of ice II is enhanced by the inclusion of atomic quantum delocalisation effects. 

In this work a re-parameterisation of the TIP4P/2005 model is provided that 'com- 
pensates' for the quantum effects so as to maintain the quantitative performance of the 
TIP4P/2005 model, whilst at the same time reproducing the correct physics at low temper- 
atures. In this new model, which we have called TIP4PQ/2005, the only parameter to have 
changed is that of the dipole moment; the charge on the hydrogen atom has been increased 
by 0.02e, thus the dipole moment increases from 2. SOD to 2.38D. 

In this paper it is shown that the TIP4PQ/2005 model provides a good description 
of the densities of the ice phases for the state points considered. The ice Ih p = I bar 
isobar has been calculated and the tendency for a to become zero is now present in the 
equation of state. This new model also correctly describes the relative stabilities of ices Ih 
and II. An extrapolation indicates that at OK Ih is more stable than ice II by 0.04 kcal/mol 
(compared to 0.014 kcal/mol experimentally). The inclusion of quantum effects substantially 
improves the overall description of all of the ice phases studied here. The TIP4P/2005 does 
a reasonable job, but the TIP4PQ/2005 is clearly superior. This paper can be regarded as 
a first step in introducing atomic quantum delocalisation effects in the description of the 
solid phases of water. However, it is by no means the last word, since obviously water is a 
flexible molecule. In our opinion the results in the present manuscript could be very useful 
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as a point of departure for the development of a flexible model of water for use in path 
integral simulations, and provides valuable material from which to make comparisons. Such 
comparison would establish just how much of the quantum effects in water are due to intra 
and how much is due to the intermolecular degrees of freedom. 
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TABLE I: Parameters for both the TIP4P/2005 and the TIP4PQ/2005 models. The distance 
between the oxygen and hydrogen sites is don- The angle, in degrees, formed by hydrogen, oxygen, 
and the other hydrogen atom is denoted by ZH-O-H. The Lennard-Jones site is located on the 
oxygen with parameters a and e. The charge on the proton is qu- The negative charge is placed 
in a point M at a distance doM from the oxygen along the H-O-H bisector. 



Model don (A) 


ZH-O-H 


a(A) 


e/kB{K) 


9H(e) 


c?om(A) 


TIP4P/2005 0.9572 


104.52 


3.1589 


93.2 


0.5564 


0.1546 


TIP4PQ/2005 0.9572 


104.52 


3.1589 


93.2 


0.5764 


0.1546 
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TABLE II: Results for the TIP4P/2005(pi) model for the systems studied, along with a comparison to classical results for the same model. 
All energies are in units of kcal/mol and the densities are in g-cm~^. The errors (in kcal/mol) are 0(0.003) in i^translationab 0(0.02) in 
^rotational, 0(0.02) in U, 0(0.04) in E and 0(0.002) g-cm'^ in p . 



Phase (& N° 


molecules) T (K) p (bars 


) (3/2) i?r -f^translational 


-^rotational 


-f^total 


U 


E 


U (classical) p (path-intej 


^ral) p (classical) 


h (432) 


250 





0.75 


0.83 


1.39 


2.22 


-12.38 


-10.17 


-13.35 


0.899 


0.920 


Ic (216) 


78 





0.23 


0.45 


1.36 


1.81 


-13.03 


-11.22 


-14.58 


0.906 


0.943 


II (432) 


123 





0.37 


0.51 


1.26 


1.77 


-12.83 


-11.06 


-14.07 


1.160 


1.198 


III (324) 


250 


2800 


0.75 


0.83 


1.35 


2.18 


-12.15 


-9.96 


-13.06 


1.141 


1.159 


IV (432) 


110 





0.33 


0.49 


1.25 


1.74 


-12.44 


-10.70 


-13.74 


1.248 


1.292 


V (504) 


237.65 


5300 


0.71 


0.80 


1.35 


2.14 


-12.19 


-10.04 


-13.21 


1.240 


1.271 


VI (360) 


225 


11000 


0.67 


0.78 


1.34 


2.12 


-12.21 


-10.10 


-13.11 


1.356 


1.379 


VII (432) 


300 


100000 


0.89 


1.05 


1.44 


2.49 


-9.32 


-6.83 


-9.95 


1.767 


1.782 


VIII (600) 


77 


24000 


0.23 


0.49 


1.17 


1.76 


-11.31 


-9.65 


-12.50 


1.573 


1.616 


IX (324) 


165 


2800 


0.49 


0.63 


1.33 


1.96 


-12.80 


-10.84 


-13.95 


1.160 


1.190 


XI (360) 


77 





0.23 


0.45 


1.36 


1.81 


-13.04 


-11.23 


-14.60 


0.906 


0.945 


XII (540) 


260 


5000 


0.77 


0.86 


1.34 


2.20 


-11.97 


-9.77 


-12.85 


1.267 


1.296 


XIII (504) 


80 


1 


0.24 


0.44 


1.25 


1.69 


-12.76 


-11.07 


-14.16 


1.217 


1.261 


XIV (540) 


80 


1 


0.24 


0.44 


1.27 


1.71 


-12.80 


-11.09 


-14.25 


1.280 


1.331 



TABLE III: Oxygen-oxygen radial distribution function of ice for various water models at 250K 
and p=0 bar. 



Model 


peak 1 


peak 2 


Reference 


r height 


r height 




TIP4P (classical) 


2.725A 4.707 


4.525A 2.279 


68 


TIP4P (path integral) 


2.7625A 4.167 


4.5625A 2.122 


This work 


TIP4P/2005 (classical) 


2.7375A 5.113 


4.5125A 2.382 


This work 


TIP4P/2005(pi) 


2.7875A 4.481 


4.5875A 2.270 


This work 


TIP4PQ/2005 


2.7625A 4.725 


4.5375A 2.405 


This work 
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TABLE IV: Results for the TIP4P/2005(pj-) model for the low temperature ice phases at a pressure 
of 1 bar. The energies are in units of kcal/mol and the densities are in g-cm~^. The errors (in 
kcal/mol) are 0(0.003) in i^transiationab 0(0.02) in i^rotationai, 0(0.02) in U, 0(0.04) in E and 
0(0.002) g-cm-=^ in p . 



Phase T (K) i^translational 


-^rotational 


total 


u 


E o 






n 7n 

VJ. 1 u 


1.36 


2.06 


-12.62 


-10 56 903 




-L(JU 


XJ.OO 


1.35 


1.93 


-12.84 


-10 91 906 




125 


0.53 


1.35 


1.87 


-12.92 


-11 05 907 




100 


0.48 


1.35 


1.83 


-12.99 


-11.15 0.907 




77 


0.45 


1.36 


1.80 


-13.02 


-11.22 0.906 




200 


n fiQ 


1.28 


1.96 


-12.54 


-10 57 1 145 




-L(JU 


XJ.O i 


1.25 


1.82 


-12.77 


-10 95 1 155 




125 


0.51 


1.24 


1.75 


-12.85 


-11 09 1 159 




100 


0.47 


1.24 


1.71 


-12.92 


-11.21 1.163 




77 


0.43 


1.26 


1.84 


-12.94 


-11.26 1.165 








1.32 


2.02 


-12.35 


-10 34 1 106 




-L(JU 


yj.oo 


1.29 


1.87 


-12.58 


-10 70 1 116 




125 


0.53 


1 30 


1.82 


-12 66 


-10 84 1 122 




mo 

J.\J\J 


n 48 


1.30 


1.77 


-12.74 


-10.96 1.125 


III 


77 


0.44 


1.30 


1.74 


-12.77 


-11.02 1.130 


V 

V 


200 


fiQ 


1.30 


1.99 


-12.28 


-10 29 1 204 


V 

V 


1 ,^0 


57 


1.28 


1.85 


-12.51 


-10 67 1 217 


V 


125 


0.52 


1.28 


1.78 


-12.60 


-10 81 1 222 


V 


100 


0.47 


1.28 


1.74 


-12.67 


-10.92 1.225 


V 


77 


0.44 


1.28 


1.72 


-12.70 


-10.99 1.227 


VI 


200 


0.69 


1.27 


1.96 


-12.19 


-10.22 1.282 


VI 


150 


0.58 


1.25 


1.83 


-12.41 


-10.58 1.296 


VI 


125 


0.51 


1.25 


1.76 


-12.50 


-10.74 1.302 


VI 


100 


0.46 


1.25 


1.71 


-12.57 


-10.86 1.306 


VI 


77 


0.43 


1.26 


1.69 


-12.60 


-10.91 1.309 
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TABLE V: Comparison of the energies, E, at OK for a selection of phases for both the 
TIP4P/2005(pi) and the TIP4PQ/2005 models as weh as results for the classical TIP4P/2005 
models. The energies are in units of kcal/mol. The lowest energy (most stable phase) is shown in 
bold font. The lower section provides the relative energies with respect to ice II. 
Ice E (OK estimate) 



TIP4P/2005 TIP4P/2005(pi) TIP4PQ/2005 Experimental^ 



h 


-15.059 


-11.240 


-12.477 


-11.315 


II 


-14.847 


-11.290 


-12.436 


-11.301 


III 


-14.741 


-11.048 


-12.210 


-11.100 


V 


-14.644 


-11.013 


-12.152 


-11.088 


VI 


-14.513 


-10.939 


-12.033 


-10.928 


Ih 


-0.212 


0.050 


-0.041 


-0.014 


II 














III 


0.106 


0.242 


0.226 


0.201 


V 


0.203 


0.277 


0.285 


0.213 


VI 


0.334 


0.351 


0.403 


0.373 
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TABLE VI: PIMC results for the TIP4PQ/2005 model for the systems studied and their relation to the experimental densities. All energies 
are in units of kcal/mol and the densities are in g-cm~'^. The errors (in kcal/mol) are 0(0.003) in i^transiationah 0(0.02) in Krotationah 0(0.02) 
in [/, 0(0.04) in E and 0(0.002) g-cm'^ in p . 



Phase T (K) p (bars 


) (3/2) i?r i^translational 


-^rotational 


-^total 


U 


E 


p (path-integral) p (experimental) Reference 




250 





0.75 


0.83 


1.45 


2.28 


-13.74 


-11.46 


0.921 


0.920 


85 


Ic 


78 





0.23 


0.46 


1.43 


1.89 


-14.33 


-12.44 


0.925 


0.931 


86 


II 


123 





0.37 


0.52 


1.32 


1.85 


-14.06 


-12.21 


1.185 


1.190 


87 


III 


250 


2800 


0.75 


0.84 


1.41 


2.25 


-13.44 


-11.18 


1.159 


1.165 


88 


IV 


110 





0.33 


0.50 


1.32 


1.82 


-13.63 


-11.81 


1.276 


1.272 


89 


V 


237.65 


5300 


0.71 


0.81 


1.41 


2.22 


-13.43 


-11.21 


1.266 


1.271 


90 


VI 


225 


11000 


0.67 


0.79 


1.39 


2.18 


-13.41 


-11.23 


1.377 


1.373 


91 


VII 


300 


100000 


0.89 


1.05 


1.47 


2.52 


-10.37 


-7.85 


1.780 


1.880 


92 


VIII 


77 


24000 


0.23 


0.50 


1.23 


1.73 


-12.28 


-10.56 


1.592 


1.628 (at lOK) 


91 


IX 


165 


2800 


0.49 


0.64 


1.39 


2.04 


-14.07 


-12.03 


1.182 


1.194 


88 


XI 


77 





0.23 


0.46 


1.43 


1.89 


-14.34 


-12.46 


0.926 


0.934 (at 5K) 


93 


XII 


260 


5000 


0.77 


0.87 


1.40 


2.27 


-13.23 


-10.96 


1.297 


1.292 


94 


XIII 


80 


1 


0.24 


0.46 


1.32 


1.77 


-13.95 


-12.17 


1.242 


1.244 


95 


XIV 


80 


1 


0.24 


0.46 


1.34 


1.80 


-13.99 


-12.20 


1.307 


1.332 


95 



TABLE VII: Unit cell parameters for the TIP4PQ/2005 model for a selection of ice phases. Ex- 
perimental values are from Table 11.2 of Ref 3. Note that for ice II the hexagonal unit cell, rather 
than the rhombohedral unit cell, is given. All distances are in angstroms. 



Phase T (K) p(bars) unit cell 



experimental simulation 





250 





a=4.518, c=7.356 


a=4.483, c= 7.352 


II 


123 





a=12.97, c=6.25 


a= 12.98, c=6.23 


III 


250 


2800 


a=6.666, c=6.936 


a=6.645, c=7.011 


V 


100 


1 


a=9.22, b =7.54, 


a=9.06, b=7.64, 








c=10.35, (3 = 109.2° 


c=10.21, p = 108.6° 


VI 


225 


11000 


a=6.181, c=5.698 


a=6.167, c=5.713 
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TABLE VIII: PIMC results for the TIP4PQ/2005 model for the low temperature ice phases at a 
pressure of 1 bar. All energies are in units of kcal/mol and the densities are in g-cm~^. The errors 
(in kcal/mol) are 0(0.003) in i^transiationai, 0(0.02) in ^rotational, 0(0.02) in U, 0(0.04) in E and 
0(0.002) g-cm-=^ in p . 



Phase T (K) i^translational 


-^rotational 


total 


u 


E o 




OuU 


n Q7 


1.49 


2.47 


-13.46 


-10 99 915 




200 


0.71 


1.44 


2.15 


-13.98 


-11 82 925 




150 


0.60 


1.42 


2.02 


-14.18 


-12.16 0.928 




125 


0.54 


1.41 


1.96 


-14.25 


-12.29 0.928 




100 


50 


1.42 


1.92 


-14.32 


-12.40 0.928 




77 


0.46 


1.43 


1.89 


-14.34 


-12 45 927 


II 


125 


0.53 


1.32 


1.84 


-14.06 


-12.21 1.185 


II 


mo 

J.\J\J 


n 48 


1.30 


1.78 


-14.14 


-12.35 1.188 


II 


77 


0.44 


1.32 


1.76 


-14.16 


-12.40 1.190 


III 


150 


5Q 


1.36 


1.95 


-13.83 


-11.88 1.134 


III 


125 


0.54 


1.36 


1.90 


-13.92 


-12.02 1.139 


III 


100 


0.50 


1.37 


1.87 


-13.98 


-12.12 1.142 


III 


77 


0.45 


1.37 


1.82 


-14.01 


-12.19 1.146 


V 


125 


0.53 


1.34 


1.84 


-13.80 


-11.93 1.248 


V 


100 


0.49 


1.34 


1.82 


-13.88 


-12.06 1.251 


V 


77 


0.44 


1.35 


1.79 


-13.91 


-12.12 1.253 


VI 


125 


0.53 


1.32 


1.85 


-13.67 


-11.82 1.330 


VI 


100 


0.48 


1.31 


1.79 


-13.74 


-11.95 1.334 


VI 


77 


0.45 


1.32 


1.77 


-13.77 


-12.00 1.336 
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TABLE IX: Estimates of the coexistence pressures (in bar) for the TIP4PQ/2005 model extrap- 
olated to OK. Experimental values are taken from the work of Whallej^ and the values for the 
classical TIP4P/2005 model are from^!^. 



Phases TIP4P/2005 TIP4PQ/2005 Experimental value 





2090 


400 


140 ± 200 


ih-ni 


3630 


3008 


2400 ± 100 


II-V 


11230 


15630 


18500 ± 4000 


II-VI 


8530 


10190 


10500 ± 1000 


III-V 


3060 


1800 


3000 ± 100 


V-VI 


6210 


5580 


6200 ± 200 
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20 40 60 80 100 120 140 

Temperature (K) 

FIG. 1: Kinetic rotational energy from PIMC simulations of the isolated H2O molecule (filled 
circles) as a function of temperature. Between 10 and 50 replicas (P) have been used, depending 
on the temperature. There is good agreement between the simulation data and the rotational 
energy obtained from the theoretical partition function of an asymmetric top having the H2O 
geometry (solid line) . The magnitude of the error is less than the size of the symbols shown. 
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2 3 4 5 6 7 

ro-0 (A) 

FIG. 2: Radial distribution function of ice Ih for TIP4P/2005 (dashed green line) and 
TIP4P/2005(pi) (solid red line) at 250 K and p=0 bar. The blue dotted line corresponds to 
the experimental data of Soper at 220K— . 
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2 3 4 5 6 7 
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FIG. 3: Radial distribution function of ice II for TIP4P/2005 (dashed green line) and 
TIP4P/2005(pi) (solid red line) at 123 K and p=0 bar. 
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FIG. 4: Radial distribution function of ice VI for TIP4P/2005 (dashed green line) and 
TIP4P/2005(pi) (solid red line) at 225 K and p=ll kbar. 
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FIG. 5: Equations of state for ice I/i at p = 1 bar. Classical TIP4P/2005 model (grey dot-dashed 
line / filled triangles)^, experimental data (red solid line)^^, TIP4P/2005(pi) (blue dotted line/ 
filled squares) and the new TIP4PQ/2005 model (black double-dotted line / filled circles). The 
error in the density is of order ±0.002 g-cm~'^. 
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20 40 60 80 100 120 

Temperature (K) 



FIG. 6: Plot of the total energy of ices Ih, II, HI, V and VI at low temperatures for p=l bar for 
TIP4P/2005(pi). Lines correspond to the fit E = a + bT^. The error in the total energy is of order 
±0.04 kcal/mol. 
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FIG. 7: Radial distribution function of ice Ih for the TIP4PQ/2005 model using PIMC (dashed 
blue line) compared with the classical TIP4P/2005 model (dotted red line) and with experimental 
data (solid red line)^ at 77K and p=l bar. 
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FIG. 8: Plot of the total energy of ices I^j, II, III, V and VI at low temperatures for p=l bar for 
the TIP4PQ/2005 model. Lines correspond to the fit = a + hT^. The error in the total energy 
is of order ±0.04 kcal/mol. 
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